Overlay Trait Data

Libraries

library(tidyverse)
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
library(readxl)

Make species x trait matrix

Greenhouse traits

• height, leaf dry matter content, specific leaf area, root mass fraction, relative growth rate, coarse root diameter, root density, specific root length of fine roots, specific root length of coarse roots

Field traits

• height, leaf dry matter content, specific leaf area, seed mass, leaf C, leaf N

“H1: Pairing field results with a greenhouse 15 N experiment will allow us to relate species responses to their ability to uptake organic N”

“In addition, we will conduct a greenhouse experiment to isolate the mechanism (mineralized N vs organic N uptake rates ) by which compost may affect species composition”

Workflow plan:

• overlay species trait data as vectors on the NMDS

  • only greenhouse traits available (only sla data from field)

    • use PC scores as composite trait axes OR

    • average the traits for each species

Read in and out matrices from traitOrdination code

Read in trait data

# # read in greenhouse trait data
# traits.gh <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=ya1ck6p7&dl=1")
# 
# # read in field trait data
# # download.file("https://www.dropbox.com/scl/fi/trqrnz54dkjoa9do4tle5/leaf-area-data_NG.xlsx?rlkey=s0kic2a5yqknikc4icek8y7rn&st=dl4x634h&dl=1",destfile = "field_sla.xlsx")
# traits.sla.field <- 
#   read_excel("field_sla.xlsx",col_names = TRUE,col_types = "text")
# 
# rda(traits.gh, scale = TRUE)

Investigate trait data, looking in Greenhouse trait folder

Trait data is identical. Will clean data in

# dropbox.main <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=0h1bs53j&dl=1")
# 
# dropbox.gh <- read.csv("https://www.dropbox.com/scl/fi/67aaubdydiy4ixffmwvfe/GreenhouseTraits.csv?rlkey=rm2eqeqome7hgoedhnal0gx6d&st=vkym6evx&dl=1")
# 
# all.equal(dropbox.main, dropbox.gh)
# identical(dropbox.main, dropbox.gh)

traits.cleaned <- read.csv("data/GreenhouseTraits_corrected_20240212.csv") %>%
  mutate(
    Fresh.leaf.mass..g. = as.numeric(ifelse(
      Fresh.leaf.mass..g. == "Skipped accidentally", NA, Fresh.leaf.mass..g.
    )),
    Dry.leaf.mass..g. = as.numeric(ifelse(
      Dry.leaf.mass..g. == "No sample", NA, Dry.leaf.mass..g.
    )),
    Root.dry.biomass..g. = as.numeric(ifelse(
      Root.dry.biomass..g. == "roots lost", NA, Root.dry.biomass..g.
    )),
    Root.volume..cm3. = as.numeric(ifelse(
      Root.volume..cm3. == "not scanned", NA, Root.volume..cm3.
    ))
  ) 

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `Root.dry.biomass..g. = as.numeric(ifelse(Root.dry.biomass..g. == "roots lost", NA, Root.dry.biomass..g.))`.
Caused by warning:
! NAs introduced by coercion

Summarize trait data (GH) for each species

means

# make trait species ID match matrix data (notes shown below) -------
traits.cleaned <- traits.cleaned %>%
  mutate(
    ID = ifelse(
      ID == "Agoseris", "AGSP", ifelse(
        ID == "AVEBAR", "AVBA", ifelse(
          ID == "BRDR", "BRDI", ID
        )
      )
    )
  )

# summary(traits.cleaned)
traits.means <- traits.cleaned %>%
  group_by(ID) %>%
  summarise(
    across(4:(ncol(traits.gh) - 1), list(mean = ~mean(.x, na.rm = TRUE)
  ))) %>%
  rename(
    Height = Height..cm._mean,
    FreshLeafMass = Fresh.leaf.mass..g._mean,
    DryLeafMass = Dry.leaf.mass..g._mean,
    LDMC = LDMC_mean, # Leaf dry matter content (LDMC, the ratio of leaf dry mass to fresh mass)
    LeafArea = Leaf.Area..cm2._mean,
    SLA = SLA..cm2.g._mean,
    ShootDryBiomass = Shoot.dry.biomass..g._mean,
    RootDryBiomass = Root.dry.biomass..g._mean,
    TotalBiomass = Total.biomass..g._mean,
    RMF = RMF_mean, # Root mass fraction (RMF) is a plant trait that measures the proportion of a plant's dry mass that is in its roots
    RootVol = Root.volume..cm3._mean,
    RootDensity = Root.density..g.cm3._mean,
    CoarseRootDiameter = Coarse.root.diameter..mm._mean,
    RootLength = Length..mm._mean,
    FineRootLength = Fine.root.length..mm._mean,
    CoarseRootLength = Coarse.root.length..mm._mean,
    CoarseRootSpecLength = Coarse.root.specific.length..cm.g._mean,
    FineRootSpecLength = Fine.root.specific.length..cm.g._mean,
    PropFineRoots = Proportion.fine.roots_mean
  )
  

# make site by species matrix match traits species
matrix.names <- data.frame(colnames(matrix.bray)) #species comp (75 count)
trait.names <- data.frame(traits.means$ID) #trait species (80 count)

intersect.names <- intersect(colnames(matrix.bray), traits.means$ID) # same species (43 count)
setdiff(colnames(matrix.bray), traits.means$ID)

 [1] "ASSP"  "AST1"  "AST2"  "AST3"  "CIQU"  "CRTI"  "CYEC"  "DACA"  "DIMU"  "GAGR"  "GAPA"  "GAPH"  "HYSP"  "JUBU"  "LYHY"  "MAD1"  "MAD2"  "NAPU"  "PAVI"  "SABI1"
[21] "SABI2" "SIGA"  "THGR"  "TRSP"  "UNBU"  "UNF1"  "UNF3"  "UNF4"  "UNF8"  "UNFR"  "UNGR"  "XAST" 

setdiff(traits.means$ID, colnames(matrix.bray))

 [1] "ACMI"      "ANARp"     "AVFA"      "BRCA"      "BRHOp"     "BRNIf"     "BRNIp"     "CLPUf"     "CLPUp"     "CYDA"      "Crassula"  "ESCA"      "FEMI"     
[14] "FEMY"      "GITRf"     "GITRp"     "HYGL"      "LACA"      "LENIf"     "LENIp"     "LOPU"      "LUBI"      "Leontodon" "MAELf"     "MAELp"     "MICA"     
[27] "NEMA"      "PLERf"     "PLERp"     "RUPU"      "Sanicula"  "TRCI"      "TRHIpi"    "TRWIf"     "TRWIp"     "TRWIpi"    "VIVA"     

# check matrix codes to see if same species are named different things
unique.matrix.species <- cover.dat %>%
  dplyr::select(code4, species) %>%
  summarise(
    code4 = unique(code4),
    species = unique(species)
  )

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.

unique.trait.species <- traits.cleaned %>%
  dplyr::select(Taxon, ID) %>%
  distinct(ID, .keep_all = TRUE)
  

# write traits.means csv ----
# write.csv(traits.means, "data/traits/all_sp_trait_means.csv")

See if I have to match names with traits ever year (ANS: All Same Species)

# # all years -------------
# # make site by species matrix match traits species
# matrix.names <- data.frame(colnames(matrix.bray)) #species comp (75 count)
# trait.names <- data.frame(traits.means$ID) #trait species (80 count)
# 
# intersect.names <- intersect(colnames(matrix.bray), traits.means$ID) # same species (43 count)
# setdiff(colnames(matrix.bray), traits.means$ID)
# setdiff(traits.means$ID, colnames(matrix.bray))
# 
# # check matrix codes to see if same species are named different things
# unique.matrix.species <- cover.dat %>%
#   dplyr::select(code4, species) %>%
#   summarise(
#     code4 = unique(code4),
#     species = unique(species)
#   )
# 
# unique.trait.species <- traits.cleaned %>%
#   dplyr::select(Taxon, ID) %>%
#   distinct(ID, .keep_all = TRUE)
#   
# # 2019 -------------
# data.frame(colnames(matrix.bray.2019)) #species comp (75 count)
# data.frame(colnames(matrix.bray.2020)) #species comp (75 count)
# data.frame(colnames(matrix.bray.2021)) #species comp (75 count)

• change Trait name: ‘Agoseris’ –> ‘Agoseris unknown’; assuming the same as ‘AGSP’ –> ‘Agoseris sp.’ in matrix data (change to matrix name)

• What’s the difference between:

  • ‘ANAR’ and ‘ANARp’ in trait data?

  • ‘BRHO’ and ‘BRHOp’ in trait data?

  • ‘BRNIf’ and ‘BRNIp’ in trait data? (not in matrix data)

  • ‘CLPUf’ and ‘CLPUp’ in trait data? (not in matrix data)

  • ‘GITRf’ and ‘GITRp’ in trait data? (not in matrix data)

  • ‘LENIf’ and ‘LENIp’ in trait data? (not in matrix data)

  • ‘TRHI’ and ‘TRHIpi’ (same sp name) in trait data? (TRHI in matrix data)

  • ‘TRWIf’ and ‘TRWIp’ and ‘TRWIpi’ in trait data? (not in matrix data)

  • ‘Crassula’ and ‘Crassula tillaea’ (could they be the same?)

  • ‘Sanicula’ (unknown) and ‘SABI1’ or ‘SABI2’ (could they be the same?)

  • ‘MAELf’ (trait data) and ‘MAD1’ (Madia sp. 1 (“tall tarweed”)’ (could they be the same?)

  • ‘MAELp’ (trait data) and ‘MAD2’ (Madia sp. 1 (“small tarweed”)’ (could they be the same?)

  • ‘Trifolium eriocephalum’ (TRER) and ‘Triphysaria eriantha’ (TRER) (could they be the same?)

• change ‘AVEBAR’ (trait data) to ‘AVBA’ (matrix data)

• change BRDR –> BRDI

Overlay Trait Vectors onto NMDS

Aboveground Traits

• Height

• Leaf

  • FreshLeafMass

  • DryLeafMass

  • LDMC

    • Leaf dry matter content (LDMC, the ratio of leaf dry mass to fresh mass)

  • LeafArea

  • SLA

  • ShootDryBiomass

Belowground Traits

• RootDryBiomass

• RMF

  • Root mass fraction (RMF) is a plant trait that measures the proportion of a plant’s dry mass that is in its roots

• RootVol

• RootDensity

• RootLength

• PropFineRoots

• Coarse Roots

  • CoarseRootDiameter

  • CoarseRootLength

  • CoarseRootSpecLength

• Fine Roots

  • FineRootLength

  • FineRootSpecLength

• TotalBiomass

ALL YEARS – LOW AND HIGH GRAZING

Prepare NMDS and Traits

Stress:

0.2011127

# subset nmds to match species names
nmds.comp.trait <- metaMDS(matrix.bray %>%
                             dplyr::select(intersect.names), 
                           k=2, trymax = 25)

Run 0 stress 0.2060005 
Run 1 stress 0.2019906 
... New best solution
... Procrustes: rmse 0.07567441  max resid 0.2449078 
Run 2 stress 0.2007527 
... New best solution
... Procrustes: rmse 0.08985478  max resid 0.2109131 
Run 3 stress 0.2094711 
Run 4 stress 0.2107847 
Run 5 stress 0.2089838 
Run 6 stress 0.2078775 
Run 7 stress 0.2027916 
Run 8 stress 0.2173911 
Run 9 stress 0.2059157 
Run 10 stress 0.2220001 
Run 11 stress 0.2079152 
Run 12 stress 0.2124366 
Run 13 stress 0.2050097 
Run 14 stress 0.202904 
Run 15 stress 0.2027187 
Run 16 stress 0.2032121 
Run 17 stress 0.2179617 
Run 18 stress 0.2028023 
Run 19 stress 0.2012207 
... Procrustes: rmse 0.02040594  max resid 0.1298012 
Run 20 stress 0.2051096 
Run 21 stress 0.2080723 
Run 22 stress 0.201363 
Run 23 stress 0.2029037 
Run 24 stress 0.2173911 
Run 25 stress 0.2079715 
*** Best solution was not repeated -- monoMDS stopping criteria:
    22: stress ratio > sratmax
     3: scale factor of the gradient < sfgrmin

nmds.comp.trait # pretty not great stress

Call:
metaMDS(comm = matrix.bray %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.2007527 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 2 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray %>% dplyr::select(intersect.names)’ 

stressplot(nmds.comp.trait)

plot(nmds.comp.trait, type="t")

nmds1 <- as.data.frame(scores(nmds.comp.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.comp.trait, choices=c(2), display=c("sites")))

nmds_dat_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# subset matrix and trait data to match names
traits.means.comp <- traits.means %>%
  filter(ID %in% intersect.names) 

# traits.means.comp2 <- traits.means.comp %>%
# dplyr::select(-ID)
rownames(traits.means.comp2) <- traits.means.comp$ID

Warning: Setting row names on a tibble is deprecated.

# trait vectors onto the NMDS
# envfit(nmds.comp.trait, traits.means.comp, permutations = 999)
# example from stack overflow: envfit(ord ~ var1 + var2 + var3, dune.spec, display="sp"), 'ord' is the nmds model, 'var1' is the trait1 i think?

colnames(traits.means.comp)

 [1] "ID"                   "Height"               "FreshLeafMass"        "DryLeafMass"          "LDMC"                 "LeafArea"             "SLA"                 
 [8] "ShootDryBiomass"      "RootDryBiomass"       "TotalBiomass"         "RMF"                  "RootVol"              "RootDensity"          "CoarseRootDiameter"  
[15] "RootLength"           "FineRootLength"       "CoarseRootLength"     "CoarseRootSpecLength" "FineRootSpecLength"   "PropFineRoots"       

vectors.nmds.comp.trait <- envfit(nmds.comp.trait ~ Height + 
                                    LeafArea +
                                    RMF, traits.means.comp, display = "sp") # ONLY Height, Leaf area, and RMF significant

#TEMPLATE: 
  # envfit(nmds.block.trait ~ Height + 
  #          FreshLeafMass +
  #          DryLeafMass +
  #          LDMC +
  #          LeafArea +
  #          SLA +
  #          ShootDryBiomass +
  #          RootDryBiomass +
  #          TotalBiomass +
  #          RMF +
  #          RootVol +
  #          RootDensity +
  #          CoarseRootDiameter +
  #          RootLength +
  #          FineRootLength +
  #          CoarseRootLength +
  #          CoarseRootSpecLength +
  #          FineRootSpecLength +
  #          PropFineRoots, traits.means.comp, display = "sp", na.rm = TRUE)

# is it different when you only do certain traits (above- vs. below-ground)? or can you just throw everything in there and it will tell you what is significant and you can subset from there?

# Looks like different subsets of traits put in are the same, so I shall proceed with the mega models and verify with Lauren

trait_vectors <- as.data.frame(vectors.nmds.comp.trait$vectors$arrows)

NMDS model: only Height, Leaf area, and RMF () found significant.

Plot

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment") +
  custom_colors +
  scale_color_manual(
    values = c(
      "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
      "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
      "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
      "wet" = rgb(0, 0, 255, maxColorValue = 255),
      "dry" = rgb(183, 65, 14, maxColorValue = 255)
    )
  ) +
  geom_segment(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) +
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year")+
  geom_segment(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) +
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)+
  geom_segment(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) +
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Grazing History") +
  geom_segment(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

# ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  color = as.factor(grazing_hist), fill = as.factor(yr))) +
#   geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
#   stat_ellipse(aes(color = as.factor(yr)),
#                 geom = "polygon",
#                fill = NA,
#                size = 1) +
#   labs(title = "Grouped by Grazing History & Year") +
#   scale_color_manual(
#   values = c(
#     "low" = "cyan",
#     "high" = "brown",
#   "2019" = "salmon",
#   "2020" = "limegreen",
#   "2021" = "cornflowerblue")
#   ) +
#   geom_segment(
#     data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
#     aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
#     inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
#     arrow = arrow(length = unit(0.2, "cm")), 
#     color = "red", 
#     size = 1
#   ) +
#   geom_text(
#     data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
#     aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)), 
#     inherit.aes = FALSE,  # Avoids unwanted inheritance
#     hjust = -0.2, 
#     vjust = -0.2, 
#     color = "red",
#     size = 5
#   )
#   
# ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
#   geom_point(shape = 21, size = 4, stroke = 1) +
#   stat_ellipse(aes(color = as.factor(nut_trt)),
#              geom = "polygon",
#              fill = NA,  # Makes the ellipse transparent inside
#              size = 1) +
#   labs(title = "Grouped by Nutrient Treatment & Year") +
#   scale_color_manual(
#   values = c(
#     "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
#     "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
#     "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
#     "wet" = rgb(0, 0, 255, maxColorValue = 255),
#     "dry" = rgb(183, 65, 14, maxColorValue = 255)
#   )
# ) +
#   geom_segment(
#     data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
#     aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
#     inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
#     arrow = arrow(length = unit(0.2, "cm")), 
#     color = "red", 
#     size = 1
#   ) +
#   geom_text(
#     data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
#     aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)), 
#     inherit.aes = FALSE,  # Avoids unwanted inheritance
#     hjust = -0.2, 
#     vjust = -0.2, 
#     color = "red",
#     size = 5
#   )

ALL YEARS – LOW AND HIGH GRAZING (Blocks)

Stress:

0.2353937

Prepare NMDS and Traits

# subset nmds to match species names
nmds.block.trait <- metaMDS(matrix.block.bray %>%
                             dplyr::select(intersect.names), , k=2, trymax = 25)

Run 0 stress 0.2353937 
Run 1 stress 0.250124 
Run 2 stress 0.2353908 
... New best solution
... Procrustes: rmse 0.002559202  max resid 0.02324147 
Run 3 stress 0.2386994 
Run 4 stress 0.2487719 
Run 5 stress 0.2453858 
Run 6 stress 0.2574031 
Run 7 stress 0.2354028 
... Procrustes: rmse 0.002945858  max resid 0.02371427 
Run 8 stress 0.2563643 
Run 9 stress 0.2460121 
Run 10 stress 0.2644327 
Run 11 stress 0.2560224 
Run 12 stress 0.246987 
Run 13 stress 0.2354028 
... Procrustes: rmse 0.002938524  max resid 0.02369694 
Run 14 stress 0.2423148 
Run 15 stress 0.2463114 
Run 16 stress 0.2354028 
... Procrustes: rmse 0.002925757  max resid 0.02366134 
Run 17 stress 0.2353936 
... Procrustes: rmse 0.002546433  max resid 0.02319281 
Run 18 stress 0.2526558 
Run 19 stress 0.2499386 
Run 20 stress 0.2480986 
Run 21 stress 0.2497376 
Run 22 stress 0.2409873 
Run 23 stress 0.2385602 
Run 24 stress 0.2521044 
Run 25 stress 0.2378706 
*** Best solution was not repeated -- monoMDS stopping criteria:
    24: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

nmds.block.trait # pretty not great stress

Call:
metaMDS(comm = matrix.block.bray %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.block.bray %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.2353908 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 2 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.block.bray %>% dplyr::select(intersect.names)’ 

nmds1 <- as.data.frame(scores(nmds.block.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.block.trait, choices=c(2), display=c("sites")))

nmds_dat_block_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# trait data
traits.means.comp 

# trait vectors onto the NMDS
vectors.nmds.block.trait <-
  envfit(nmds.block.trait ~ Height + 
           DryLeafMass +
           LeafArea +
           RootDryBiomass +
           RMF +
           RootDensity, traits.means.comp, display = "sp")
# ONLY Height, DryLeafMass, LeafArea, RMF significant (RootDryBiomass, RootDensity marginal)

trait_vectors <- as.data.frame(vectors.nmds.block.trait$vectors$arrows)

Plot

ggplot(nmds_dat_block_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment - Blocks w Traits") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_block_trait, aes(NMDS1, NMDS2,  fill = ppt_trt, color = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment - Blocks w Traits") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

NA
NA
NA
NA

All years - LOW GRAZING ONLY

Prepare NMDS and Traits

Stress:

0.1879671

# subset nmds to match species names
nmds.low.trait <- metaMDS(matrix.bray.low %>%
                             dplyr::select(intersect.names), , k=2, trymax = 25)

Run 0 stress 0.2032274 
Run 1 stress 0.1911123 
... New best solution
... Procrustes: rmse 0.07305776  max resid 0.1904385 
Run 2 stress 0.2933553 
Run 3 stress 0.1890672 
... New best solution
... Procrustes: rmse 0.1307018  max resid 0.279877 
Run 4 stress 0.2104987 
Run 5 stress 0.208207 
Run 6 stress 0.2170339 
Run 7 stress 0.2019759 
Run 8 stress 0.1911123 
Run 9 stress 0.2124695 
Run 10 stress 0.2032282 
Run 11 stress 0.2020341 
Run 12 stress 0.1890672 
... Procrustes: rmse 3.79248e-06  max resid 8.350358e-06 
... Similar to previous best
Run 13 stress 0.2138031 
Run 14 stress 0.3881609 
Run 15 stress 0.2019912 
Run 16 stress 0.1915469 
Run 17 stress 0.2112419 
Run 18 stress 0.2084351 
Run 19 stress 0.1911123 
Run 20 stress 0.2105596 
*** Best solution repeated 1 times

nmds.low.trait 

Call:
metaMDS(comm = matrix.bray.low %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.low %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.1890672 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 3 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.low %>% dplyr::select(intersect.names)’ 

nmds1 <- as.data.frame(scores(nmds.low.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.low.trait, choices=c(2), display=c("sites")))

nmds_dat_low_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# trait data: traits.means.comp 

# trait vectors onto the NMDS
vectors.nmds.low.trait <-
  envfit(nmds.low.trait ~ FreshLeafMass +
           RMF, traits.means.comp, display = "sp", na.rm = TRUE)
# ONLY FreshLeafMass (RMF marginal)

trait_vectors <- as.data.frame(vectors.nmds.low.trait$vectors$arrows)

Plot

ggplot(nmds_dat_low_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019-2021 (LOW) w Traits - Grouping by Nutrient Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  custom_colors + 
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_low_trait, aes(NMDS1, NMDS2,  fill = ppt_trt, color = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019-2021 (LOW) w Traits - Grouping by Precipitation Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  custom_colors + 
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

2019

Stress:

0.06290431

Prepare NMDS and Traits

# subset nmds to match species names
nmds.2019.trait <- metaMDS(matrix.bray.2019 %>%
                             dplyr::select(intersect.names), , k=2, trymax = 25)

Run 0 stress 0.06471496 
Run 1 stress 0.06931106 
Run 2 stress 0.06874971 
Run 3 stress 0.3616162 
Run 4 stress 0.06264873 
... New best solution
... Procrustes: rmse 0.06539037  max resid 0.2367651 
Run 5 stress 0.09229892 
Run 6 stress 0.06471497 
Run 7 stress 0.06290424 
... Procrustes: rmse 0.005564852  max resid 0.0187599 
Run 8 stress 0.06816428 
Run 9 stress 0.06901458 
Run 10 stress 0.06471496 
Run 11 stress 0.07109362 
Run 12 stress 0.06459268 
Run 13 stress 0.08080796 
Run 14 stress 0.07203292 
Run 15 stress 0.06459267 
Run 16 stress 0.07780829 
Run 17 stress 0.07518585 
Run 18 stress 0.06264873 
... Procrustes: rmse 1.714225e-05  max resid 4.819798e-05 
... Similar to previous best
Run 19 stress 0.06931106 
Run 20 stress 0.07273526 
*** Best solution repeated 1 times

nmds.2019.trait 

Call:
metaMDS(comm = matrix.bray.2019 %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2019 %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.06264873 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 4 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2019 %>% dplyr::select(intersect.names)’ 

nmds1 <- as.data.frame(scores(nmds.2019.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2019.trait, choices=c(2), display=c("sites")))

nmds_dat_2019_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# trait data: traits.means.comp 


# trait vectors onto the NMDS
vectors.nmds.2019.trait <-
  envfit(nmds.2019.trait ~ Height + 
           DryLeafMass +
           LeafArea +
           RMF +
           RootDensity +
           CoarseRootDiameter, traits.means.comp, display = "sp", na.rm = TRUE)

Plot

ggplot(nmds_dat_2019_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019 w Traits - Grouping by Nutrient Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_2019_trait, aes(NMDS1, NMDS2,  fill = ppt_trt, color = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019 w Traits - Grouping by Precipitation Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

2020

Stress:

0.1512864

Prepare NMDS and Traits

# subset nmds to match species names
nmds.2020.trait <- metaMDS(matrix.bray.2020 %>%
                             dplyr::select(intersect.names), , k=2, trymax = 25)

Run 0 stress 0.1512864 
Run 1 stress 0.16382 
Run 2 stress 0.1740134 
Run 3 stress 0.1512864 
... New best solution
... Procrustes: rmse 1.160478e-05  max resid 2.59613e-05 
... Similar to previous best
Run 4 stress 0.151653 
... Procrustes: rmse 0.01975231  max resid 0.06777473 
Run 5 stress 0.16382 
Run 6 stress 0.1601231 
Run 7 stress 0.151336 
... Procrustes: rmse 0.01055512  max resid 0.03081296 
Run 8 stress 0.1600829 
Run 9 stress 0.1922664 
Run 10 stress 0.171301 
Run 11 stress 0.1512864 
... New best solution
... Procrustes: rmse 1.626593e-05  max resid 3.63373e-05 
... Similar to previous best
Run 12 stress 0.151336 
... Procrustes: rmse 0.0105455  max resid 0.03079015 
Run 13 stress 0.1600829 
Run 14 stress 0.163364 
Run 15 stress 0.1512864 
... Procrustes: rmse 8.079558e-06  max resid 2.224799e-05 
... Similar to previous best
Run 16 stress 0.1740133 
Run 17 stress 0.1516619 
... Procrustes: rmse 0.01548936  max resid 0.05232438 
Run 18 stress 0.1714389 
Run 19 stress 0.175096 
Run 20 stress 0.16382 
*** Best solution repeated 2 times

nmds.2020.trait 

Call:
metaMDS(comm = matrix.bray.2020 %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2020 %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.1512864 
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 11 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2020 %>% dplyr::select(intersect.names)’ 

nmds1 <- as.data.frame(scores(nmds.2020.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2020.trait, choices=c(2), display=c("sites")))

nmds_dat_2020_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# trait data: traits.means.comp 


# trait vectors onto the NMDS
vectors.nmds.2020.trait <-
  envfit(nmds.2020.trait ~ FreshLeafMass +
           LeafArea, traits.means.comp, display = "sp", na.rm = TRUE)

# ONLY FreshLeafMass, LeafArea

trait_vectors <- as.data.frame(vectors.nmds.2020.trait$vectors$arrows)

Plot

ggplot(nmds_dat_2020_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 w Traits - Grouping by Nutrient Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_2020_trait, aes(NMDS1, NMDS2,  fill = ppt_trt, color = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 w Traits - Grouping by Precipitation Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

2021

Prepare NMDS and Traits

Stress:

0.123061

# subset nmds to match species names
nmds.2021.trait <- metaMDS(matrix.bray.2021 %>%
                             dplyr::select(intersect.names), , k=2, trymax = 25)

Run 0 stress 0.1234883 
Run 1 stress 0.123061 
... New best solution
... Procrustes: rmse 0.04524088  max resid 0.133887 
Run 2 stress 0.1556472 
Run 3 stress 0.1556472 
Run 4 stress 0.1510891 
Run 5 stress 0.133126 
Run 6 stress 0.1412774 
Run 7 stress 0.1234883 
... Procrustes: rmse 0.04524054  max resid 0.133303 
Run 8 stress 0.1234883 
... Procrustes: rmse 0.04524098  max resid 0.13329 
Run 9 stress 0.1331269 
Run 10 stress 0.1412774 
Run 11 stress 0.1331264 
Run 12 stress 0.1329969 
Run 13 stress 0.123061 
... Procrustes: rmse 7.588273e-06  max resid 2.28136e-05 
... Similar to previous best
Run 14 stress 0.1412774 
Run 15 stress 0.1412774 
Run 16 stress 0.1234883 
... Procrustes: rmse 0.04524105  max resid 0.1332899 
Run 17 stress 0.1234883 
... Procrustes: rmse 0.04524004  max resid 0.1333063 
Run 18 stress 0.1790971 
Run 19 stress 0.1513012 
Run 20 stress 0.123061 
... Procrustes: rmse 2.514183e-06  max resid 7.876099e-06 
... Similar to previous best
*** Best solution repeated 2 times

nmds.2021.trait 

Call:
metaMDS(comm = matrix.bray.2021 %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2021 %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.123061 
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 1 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2021 %>% dplyr::select(intersect.names)’ 

nmds1 <- as.data.frame(scores(nmds.2021.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2021.trait, choices=c(2), display=c("sites")))

nmds_dat_2021_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# trait data: traits.means.comp 
rownames(traits.means.comp) <- traits.means.comp$ID

Warning: Setting row names on a tibble is deprecated.

# trait vectors onto the NMDS
vectors.nmds.2021.trait <-
  envfit(nmds.2021.trait ~ Height + 
           FreshLeafMass +
           LeafArea, na.omit(traits.means.comp), display = "sp", na.rm = TRUE)

# (Height, FreshLeafMass, LeafArea marginal)

trait_vectors <- as.data.frame(vectors.nmds.2021.trait$vectors$arrows)

Plot

ggplot(nmds_dat_2021_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 w Traits - Grouping by Nutrient Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_2021_trait, aes(NMDS1, NMDS2,  fill = ppt_trt, color = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 w Traits - Grouping by Precipitation Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

2021 – Blocks

Stress:

0.170799

Prepare NMDS and Traits

# subset nmds to match species names
nmds.2021.block.trait <- metaMDS(matrix.bray.block.2021 %>%
                             dplyr::select(intersect.names), , k=2, trymax = 25)

Run 0 stress 0.170799 
Run 1 stress 0.2052476 
Run 2 stress 0.170799 
... New best solution
... Procrustes: rmse 9.313132e-06  max resid 3.473702e-05 
... Similar to previous best
Run 3 stress 0.2060774 
Run 4 stress 0.1873003 
Run 5 stress 0.2046024 
Run 6 stress 0.206467 
Run 7 stress 0.2180152 
Run 8 stress 0.1709819 
... Procrustes: rmse 0.007084111  max resid 0.03059497 
Run 9 stress 0.1797966 
Run 10 stress 0.170799 
... Procrustes: rmse 8.637585e-06  max resid 3.219545e-05 
... Similar to previous best
Run 11 stress 0.1797966 
Run 12 stress 0.1914994 
Run 13 stress 0.1951365 
Run 14 stress 0.1996107 
Run 15 stress 0.196785 
Run 16 stress 0.170799 
... New best solution
... Procrustes: rmse 6.097064e-06  max resid 2.279654e-05 
... Similar to previous best
Run 17 stress 0.214524 
Run 18 stress 0.1709819 
... Procrustes: rmse 0.007078194  max resid 0.03057796 
Run 19 stress 0.189269 
Run 20 stress 0.1939681 
*** Best solution repeated 1 times

nmds.2021.block.trait 

Call:
metaMDS(comm = matrix.bray.block.2021 %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.block.2021 %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.170799 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 16 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.block.2021 %>% dplyr::select(intersect.names)’ 

nmds1 <- as.data.frame(scores(nmds.2021.block.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2021.block.trait, choices=c(2), display=c("sites")))

nmds_dat_2021_block_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# trait data: traits.means.comp 


# trait vectors onto the NMDS
vectors.nmds.2021.block.trait <-
  envfit(nmds.2021.block.trait ~ DryLeafMass +
           LeafArea, traits.means.comp, display = "sp", na.rm = TRUE)

# ONLY DryLeafMass, LeafArea

trait_vectors <- as.data.frame(vectors.nmds.2021.block.trait$vectors$arrows)

Plot

ggplot(nmds_dat_2021_block_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block) w Traits - Grouping by Nutrient Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_2021_block_trait, aes(NMDS1, NMDS2,  fill = ppt_trt, color = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block) w Traits - Grouping by Precipitation Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

2021 – Blocks (LOW)

Stress:

0.1299711

Prepare NMDS and Traits

# subset nmds to match species names
nmds.2021.block.low.trait <- metaMDS(matrix.bray.block.2021.low %>%
                             dplyr::select(intersect.names), , k=2, trymax = 25)

Run 0 stress 0.1299711 
Run 1 stress 0.1299712 
... Procrustes: rmse 5.590124e-05  max resid 0.0001635787 
... Similar to previous best
Run 2 stress 0.1299712 
... Procrustes: rmse 4.398462e-05  max resid 0.0001319322 
... Similar to previous best
Run 3 stress 0.1627521 
Run 4 stress 0.1299711 
... Procrustes: rmse 1.009167e-05  max resid 3.112151e-05 
... Similar to previous best
Run 5 stress 0.162752 
Run 6 stress 0.1299711 
... Procrustes: rmse 1.406545e-05  max resid 4.211068e-05 
... Similar to previous best
Run 7 stress 0.1299711 
... Procrustes: rmse 7.661138e-06  max resid 1.674758e-05 
... Similar to previous best
Run 8 stress 0.18467 
Run 9 stress 0.2001228 
Run 10 stress 0.1602444 
Run 11 stress 0.1299712 
... Procrustes: rmse 4.253368e-05  max resid 0.0001252648 
... Similar to previous best
Run 12 stress 0.1431163 
Run 13 stress 0.2111627 
Run 14 stress 0.1431163 
Run 15 stress 0.2104642 
Run 16 stress 0.1299712 
... Procrustes: rmse 3.767893e-05  max resid 0.0001028459 
... Similar to previous best
Run 17 stress 0.1431163 
Run 18 stress 0.1785202 
Run 19 stress 0.1299711 
... New best solution
... Procrustes: rmse 8.487148e-06  max resid 2.52772e-05 
... Similar to previous best
Run 20 stress 0.1602445 
*** Best solution repeated 1 times

nmds.2021.block.low.trait 

Call:
metaMDS(comm = matrix.bray.block.2021.low %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.block.2021.low %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.1299711 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 19 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.block.2021.low %>% dplyr::select(intersect.names)’ 

nmds1 <- as.data.frame(scores(nmds.2021.block.low.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2021.block.low.trait, choices=c(2), display=c("sites")))

nmds_dat_2021_block_low_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# trait data: traits.means.comp 


# trait vectors onto the NMDS
vectors.nmds.2021.block.low.trait <-
  envfit(nmds.2021.block.low.trait ~ FreshLeafMass, traits.means.comp, display = "sp", na.rm = TRUE)

# ONLY FreshLeafMass

trait_vectors <- as.data.frame(vectors.nmds.2021.block.low.trait$vectors$arrows)

Plot

ggplot(nmds_dat_2021_block_low_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; LOW) w Trait - Grouping by Nutrient Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_2021_block_low_trait, aes(NMDS1, NMDS2,  fill = ppt_trt, color = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; LOW) w Trait - Grouping by Precipitation Treatment") + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
) +
  geom_segment(
    data = trait_vectors, 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = trait_vectors, 
    aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

2021 – Blocks (HIGH)

Prepare NMDS and Traits

# subset nmds to match species names
nmds.2021.block.high.trait <- metaMDS(matrix.bray.block.2021.high %>%
                             dplyr::select(intersect.names), , k=2, trymax = 25)

Run 0 stress 0.1416462 
Run 1 stress 0.1595099 
Run 2 stress 0.1867893 
Run 3 stress 0.2747088 
Run 4 stress 0.1867893 
Run 5 stress 0.1416462 
... New best solution
... Procrustes: rmse 3.289538e-05  max resid 8.812783e-05 
... Similar to previous best
Run 6 stress 0.1867893 
Run 7 stress 0.1416461 
... New best solution
... Procrustes: rmse 0.0004825076  max resid 0.00160539 
... Similar to previous best
Run 8 stress 0.1416461 
... Procrustes: rmse 0.0004023531  max resid 0.00133908 
... Similar to previous best
Run 9 stress 0.1595099 
Run 10 stress 0.1867893 
Run 11 stress 0.15951 
Run 12 stress 0.1867893 
Run 13 stress 0.1416461 
... Procrustes: rmse 7.442382e-05  max resid 0.000245897 
... Similar to previous best
Run 14 stress 0.1416461 
... Procrustes: rmse 4.391587e-05  max resid 0.000145225 
... Similar to previous best
Run 15 stress 0.1595099 
Run 16 stress 0.141646 
... New best solution
... Procrustes: rmse 0.0002924508  max resid 0.0009736183 
... Similar to previous best
Run 17 stress 0.1595099 
Run 18 stress 0.1416461 
... Procrustes: rmse 0.000208078  max resid 0.0006900715 
... Similar to previous best
Run 19 stress 0.141646 
... Procrustes: rmse 0.0002689142  max resid 0.0008954326 
... Similar to previous best
Run 20 stress 0.1416463 
... Procrustes: rmse 0.0002338219  max resid 0.0007719694 
... Similar to previous best
*** Best solution repeated 4 times

nmds.2021.block.high.trait 

Call:
metaMDS(comm = matrix.bray.block.2021.high %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.block.2021.high %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.141646 
Stress type 1, weak ties
Best solution was repeated 4 times in 20 tries
The best solution was from try 16 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.block.2021.high %>% dplyr::select(intersect.names)’ 

nmds1 <- as.data.frame(scores(nmds.2021.block.high.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2021.block.high.trait, choices=c(2), display=c("sites")))

nmds_dat_2021_block_high_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# trait data: traits.means.comp 


# trait vectors onto the NMDS
vectors.nmds.2021.block.high.trait <-
  envfit(nmds.2021.block.high.trait ~ Height + 
           FreshLeafMass +
           DryLeafMass +
           LDMC +
           LeafArea +
           SLA +
           ShootDryBiomass +
           RootDryBiomass +
           TotalBiomass +
           RMF +
           RootVol +
           RootDensity +
           CoarseRootDiameter +
           RootLength +
           FineRootLength +
           CoarseRootLength +
           CoarseRootSpecLength +
           FineRootSpecLength +
           PropFineRoots, traits.means.comp, display = "sp", na.rm = TRUE)

# None

trait_vectors <- as.data.frame(vectors.nmds.2021.block.high.trait$vectors$arrows)

Plot

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; HIGH) - Grouping by Nutrient Treatment") +
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  fill = ppt_trt, color = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; HIGH) - Grouping by Precipitation Treatment") +
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)